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We study how energy transport in an integrable system is affected by the spectral densities of 
heat reservoirs. The model investigated here is the quantum harmonic chain whose both ends are 
in contact with two heat reservoirs at different temperatures. The master equation for the reduced 
density matrix is derived on the assumption that the reservoirs are composed of an infinite number of 
independent harmonic oscillators. We evaluate temperature profile and energy flux in the stationary 
state for the master equation and discuss how they depend on the types of spectral densities. When 
we attach the reservoirs of the same type of spectral density, we find that the temperature profile 
is independent of the types. On the other hand, when the two reservoirs have different types of 
spectral densities, the energy profile near the ends of the chain depends on the types. When the 
coupling is finite, the temperature profile near the ends shows wide variation of behavior dependent 
on spectral densities and temperatures of reservoirs. This dependence is discussed with the Fokker- 
Planck equations obtained in the classical limit. 
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I. INTRODUCTION 



Generally integrable systems show abnormal energy transport, namely the Fourier heat law is not realized there. 
This attributes to the lack of scattering between modes, which should be induced by nonintegrability. Two typical 
characteristics are seen in the energy transport in integrable systems. One is energy flux per unit volume independent 
of the system size. The other is a flat temperature profile with no global temperature gradient. 

The harmonic chain is a typical integrable system which shows these characteristics J^-j?]]. Rieder, Lebowitz and Lieb 
(RLL) investigated the classical harmonic chain whose ends are in contact with heat reservoirs at different temperatures 
||. They exactly evaluated the covariance matrix of the variables in the stationary state using the Langevin equation, 
and they proved these characteristics. That is, they found that energy flux per volume is proportional to only 
temperature difference and is independent of the system size, and no global temperature gradient is formed. Although 
the temperature profile is flat in the internal region, they found the peculiar behavior in the vicinity of the ends of 
the chain. Namely, the local temperature is higher than the bulk value near the colder reservoir, and lower near the 
hotter reservoir. 

Ziircher and Talkner (ZT) Q] investigated a quantum model corresponding to that of RLL with use of the quantum 
Langevin equation ||. As for the bulk behavior, they found the same features as in the classical case. That is, 
no global temperature gradient is formed and energy flux is independent of system size. In the high temperature 
limit, the quantum Langevin equation is reduced to the Langevin equation with the Gaussian white noise and all the 
characteristics obtained in || are reproduced. However, temperature profile in the vicinity of the ends of the system 
shows some variety depending on temperature and a damping constant. 

The reservoir employed in these studies is only the Ohmic type, which is one of the possible three types of heat 
reservoirs: sub-Ohmic, Ohmic, and super-Ohmic. Because integrable systems have no scattering between modes, 
their nonequilibrium behavior will be easily affected by the types of reservoirs at the boundary. Thus in this paper, 
we investigate how nonequilibrium nature in the harmonic chain depends on the types of reservoirs. Here we derive 
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the master equation for the reduced density matrix through the projection operator method on assumption that the 
reservoir is composed of an infinite number of independent harmonic oscillators. This method is valid in the weak 
coupling limit and for slow motion of the system because it treats the second order perturbation with respect to a 
coupling constant and the Markovian approximation || . 

We investigate the effect of the spectral density on the temperature profile and energy flux in the quantum harmonic 
chain in contact with two reservoirs applying the master equation for reduced density matrix. At the weak coupling 
limit, the stationary state can be obtained analytically, and it is found that when the two reservoirs have the same 
type of spectral density, the temperature profile is independent of it, and the profile is the same as the result previously 
reported H[|. On the other hand, when the spectral densities are of different types, the temperature profile depends 
on the types. Even in the classical limit, the internal temperature deviates from the mean value of the temperature of 
the reservoirs and we observe a deviation of temperatures around ends from the temperature of internal region. The 
direction of deviation is determined by frequency dependence of the spectral density near uj = 

We also investigate the reduced density matrix with finite values of the coupling constant numerically. Although it 
is derived in the perturbation of the coupling constant and it is only valid in the weak coupling limit, we dare regard 
the master equation with a finite coupling constant as a model for a time evolution with a dissipation. In other words, 
we assume that the time evolution qualitatively represents a kind of real phenomena in nature, where the coupling 
constant represents the strength of the dissipative mechanism. We study qualitatively the dependence of temperature 
profile on the types of reservoirs. It even has been reported in some cases that the reduced density matrix with a 
finite coupling can produce a good long time behavior of the system in comparison with the exact path-integral result 
not only qualitatively but also quantitatively Q . 

Temperature profile at the ends of the chain is found to depend on the spectral density and the temperatures of 
the reservoirs even when reservoirs at both ends have same type of spectral density. The dependence on the types of 
the reservoirs is discussed in the Fokker-Planck equations in the classical limit. 

This paper is organized as follows. In Sec. II, we derive the master equation for the reduced density matrix of 
a general many body system in contact with a phonon reservoir. Sec. Ill is devoted to the investigation on energy 
transport in the harmonic chain coupled to the phonon reservoirs at weak coupling limit. In Sec. VI, we consider the 
case of finite coupling and investigate corresponding Fokker-Planck equations. Summary and discussions are given in 
Sec. V. 



II. THE PHONON RESERVOIR 



A. Master equation for reduced density matrix 



In this subsection, we derive the master equation for the reduced density matrix of a system in contact with a 
phonon reservoir. Let us consider the following total Hamiltonian -f/tot, 

-fftot = H + Hi nt + Hr, (2.1) 

where H denotes the Hamiltonian for the system of interest, Hn denotes the Hamiltonian for the reservoir, and Hi nt 
describes the interaction between the system and the reservoir. We assume that the reservoir consists of an infinite 
number of mutually independent harmonic oscillators |ll|-|T3]], that is, 

*=E^^=E»-(ft + i), (2.2, 

where b\ and b a are the creation and annihilation operators for the ath mode. We assume a linear coupling between 
a Hermitian operator of the system X and reservoir's operators in the form. 

a a. 

= \VhJ2l*(bl + b a )X + \' 2 Y,j^X 2 (A'>A), (2.3) 

a a 

where A is a coupling constant and 7 Q s and A' are some constants. We put the second term in the right hand side in 
order to make the total Hamiltonian to be bounded p3. This term is regarded as a part of H: 
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H 



(2.4) 



In this section, we do not make any assumption for the system, though we consider a harmonic chain for the system 
in the next section. 

We derive the master equation for the reduced density matrix following the standard method p|. We start from 
the quantum Liouville equation for the total system 



dptotjt) _ 1 ,„ M , 

57 — "TT -"totjPtot WJj 

at in 



(2.5) 



where ptot(^) is the density matrix for the total system. Under the condition that the reservoir is initially in the 
equilibrium state at inverse temperature /?, the degrees of freedom of the reservoir are traced out with the aid of 
projection operators. In order to obtain an equation which can be solved practically, we usually expand it up to 
the second order with respect to A and also adopt the Markovian approximation, which is valid when correlations 
between reservoir's variables are short-lived. As the result we obtain an equation for the reduced density matrix 
p(t) = TrR/9tot(£) (TrR means the trace concerning the reservoir's degrees of freedom) of the form 



where Tp(t) is given by 

Tp(t) = -2 dt' diue iut ' §{u){XX{-t')p{t) 
n Jo J-00 

- e 0H "Xp(t)X(-t') + e 0Ru p(t)X(-t')X - X{-t')p{t)X) . 
In Eq. (2.7), X(-t') means the Heisenberg operator at time — t' 

X(-t') = e- iHt '/ h Xe im '/ h , 



(2.6) 



(2.7) 



(2.8) 



and the function $(tu) denotes the Fourier transform of the two-point function of the reservoir's operators coupling 
to X, namely, 



where $(£) is given by 

- Tr R 



£ 



1/2 



1 

2^ 



(2.9) 



(^^) 1/2 7 Q 7 Q '^(Q)^(i)e-^ B 



Tr R e 



a 

Hence, denoting the reservoir's density of states with respect to frequency w by D(u>), we can write $(w) as 

$(w) = hj((j)^ 



. 2 D(lj)-D(-u) 



(2.10) 



(2.11) 



e l3hw _ ^ ' 

where we introduced a smooth function 7(0;) that satisfies 7(±w Q ) = 7^. Here wc define the spectral density I(to) as 

7(w) = j{lo) 2 D(lv). (2.12) 

The following form for the spectral density is considered in the literature, 

I{uj) = I a uj a 9(u;), (2.13) 

where 0(u>) is the step function: 6(u>) = lforw > and 0(u>) — Qforw < 0. The reservoir is called Ohmic if a = 1, 
sub-Ohmic if a < 1, and super-Ohmic if a > 1 [12|. 
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In the following we rewrite Eq. ( |2.6| ) in a form convenient for later use. Let us consider the matrix components of 
operator Tp(t), (k\Tp(t)\n), where \k) and \n) are eigenstates for the system Hamiltonian H with energy eigenvalues 
Ek and E n , respectively. For the integral with respect to t' , we use the mathematical formula 



e lut dt 



7n5(y) 



v 1 -, 

V 



neglecting the principal value [Q,[ij,[l8| and the Kubo-Martin-Schwinger (KMS) condition <^{u})e^ hu} 
the matrix components of operator Tp(t) is written as 



(2.14) 
Then, 



(k\r P (t)\n) = ^Yl 

l,m L 



E, 



Ei — E, 



(t) - X k ,ip hm {t)X 



E n — En 



+p k>l (t)X* l $ 



E, 



X,, 



x k ^ 



Eh — Ei 



(t)X 



Now we introduce the operator R whose matrix elements are 



(l\R\m) = \x hm § 



El — E n 



Ei-E„ 



e 0(E t -E m ) _ i 



(2.15) 



(2.16) 



Then Tp(t) is written in the following compact form, 

Tp(t) = | (XRp(t) - Rp{t)X - Xp{t)R) + p(t)R f X) 
= ~([X,Rp(t)] + [X,Rp(ttf 
Thus we arrive at the master equation of the form 

_ = _[ff,p(t)]__ 



{X,Rp(t)] + {X,Rp(ttf 



(2.17) 



(2.18) 



This is a generalized Lindblad form |15|Jl6[| treating general many body system with the coupling form (2.3). When 
the system has many body interaction, the noncommutabilities cause the operator R to contain all degrees of freedom 
of the system even if H\ nt is a part of th e system. Thus, in general R has a complicated form with all degrees freedom 
of the system. Never thel ess, Eq. ( 2.16 ) gives the explicit and compact fo rm of R for the general systems when the 
reservoir is given by (2.2). Thus we can expect that the master equation (2.18) is widely applicable for the concrete 
studies of many body systems. 

In the present study this ma ster e quation is used as a basic equation for a system coupled with the phonon reservoir. 
It is readily checked that Eq. ( [2.16 ) satisfies at least a necessary condition for the master equation, i.e., the canonical 
distribution e^^/Tr^-^ ) into p in Eq. ( |2.18| ) gives a stationary solution We also expect the stability of the 
stationary solution at least when A is small enough. 



B. Comparison with the Quantum Langevin dynamics 

Here we briefly review another type of equation representing quantum dynamics with a thermal environment that is 
called a quantum Langevin equa tion which was used in previous studies [Q and compare it with the master equation 
for reduced density matrix ( [2.18 ) (see also |l4j| for other types of quantum Langevin equations we no not explain here 

)■ 

The quantum Langevin equation was introduced by Ford, Kac, and Mazur J8|. They considered the following 
coupled oscillators composed of (27V + 1) particles, 

1 N 1 N 

H =2 ^2 p ^ + 2 51 9mA m , n qn, (2-19) 

n— — N m,n= — N 

where q n and p n are the ith canonical coordinate and momentum variable, respectively. The matrix A = {A mjl ) is a 
(2N + 1) x (2/V + 1) symmetric matrix whose elements are 
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1 N 

2N+1 ^ k 



exp 



k=-N 



. 27Tfc 

z (m — n) 

2N+V ' 



(2.20) 



It should be noted that the eigenvalues of this matrix are w 2 (s — ~N, —N + 1, • • • , N — 1, N). The authors assumed 
that the initial state of the system ( 2.19Q is in equilibrium at a temperature, and examined under what condition the 
behavior of particle can be described by a Langevin equation. 

They found the following. If the eigen- frequencies of the whole system, u s , have the special form 



f 2 tan 2 



2iV + l 



the motion of a particle of the system in the equilibrium state is described by 

dq (t) 



dt 
dppjt) 
dt 



Po, 

= -f Po + E(t), 



(2.21) 

(2.22a) 
(2.22b) 



where qo,po and E(t) are operators in the Heisenberg picture. The operator E{t) is described by operators of particles. 
In the equilibrium state, E(t) behaves like the Gaussian random force with vanishing mean (E(t)) — 0, where (...) 
means Tr (exp(— (3 H ({^(0)}, fe(0)}))...) /Z. It also satisfies the commutation relation 



[E(t),E(s)]=2ihf^- t 6(t-s) 



and has the symmetrized correlation 
1 



(E(t)E(t + t) + E(t + r)E(t)) = 



nf_ 


/>oo 

/ lu coth 


~f3hu~ 


TT 


Jo 


2 



cos(ujT)duj. 



(2.23) 



(2.24) 



This dynamics yields a classical Langevin equation with Gaussian white noise in the classical limit h — > 0. 

This dynamics has been applied to the quantum harmonic chain and investigated some quantum effects in energy 
transport phenomena ^ff). However strictly speaking this quantum Langevin equation is the dynamics for a simple 
particle system in an equilibrium state. Therefore this dynamics is not consistent with a nonequilibrium dynamics for 
many body system in principle. The master equation for reduced density matrix ( 2.1 8| ) is derived for gene ral many 
body system on the assumption that only the reservoir is in equilibrium. Thus the master equation ( 2. IS ) is more 
suitable in this context. 



III. ENERGY TRANSPORT IN THE QUANTUM HARMONIC CHAIN 

AT THE WEAK COUPLING LIMIT 



In this section, we investigate energy transport in the quantum harmonic chain in contact with two phonon reservoirs 
at different temperatures with various types of spectral density of the thermal reservoir and examine what is common 
with and what is different from the results in the classical case || and also the quantum case Q with the Ohmic 
spectral density. We first discuss the case of weak coupling limit. 



A. System 



Here we take the one-dimensional quantum harmonic chain, 

N 9 N , 

H= l^^y— + 2^ ~ x n) > (3-1) 

n—l n—0 

as the system. This Hamiltonian should be considered to be the renormalized Hamiltonian including the second term 
m (|]|. As in |], we impose the fixed boundary condition, xq = x^ + i ~ 0. By Fourier transformation: 
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X n 



2 ,. , / 2 



N , , - , - 



-^Ufcsin(fcn), p ra = W ^ 1 X! ^ sin (M> ( 3 - 2 ) 



the Hamiltonian is decoupled into the normal modes as 

majtu 



" = E^ + ^. (3-3) 



k 



where cot = 2o>o sin(/c/2). The wave number k runs through the values k = = 1) 2, • • • , N— 1, AT). It is easily 

found that operators itfc and iv satisfy the commutation relations for canonical variables 

[uk,Vk>]=iMk,k' and [u k , uy] = [v k , v k >] = 0, (3.4) 

and 

Introducing the creation and annihilation operators at and in the ordinary manner: 



mwosin(fc/2) / ivk \ + muJo sin(fc/2) / ivk 



= \> : : I u k + n . ,, /n s I and at = \/ r ' ' - I Mfc . , . . 

2mw sm(fc/2)/ * V 7i \ 2mw sm(fc/2) 

we obtain 

~ " (3-5) 



B. Equation of Motion of the System which Contacts with Two Different Reservoirs 



In order to describe the system whose ends are attached to phonon reservoirs at different temperatures, we set 
dynam ical m odel where the contacts with thermal baths are taking into account by the dissipation terms of the forms 
in Eq.(2.1S). That is, variables at the left-end and right-end points x% and xm are linearly coupled with one phonon 
reservoir at inverse temperature /3l and /3r, respectively. 

We assume that the coupling strength A and the form of the coupling 7^ in Eq. (2.3) are common for both the 
reservoirs. Then the master equation for the reduced density matrix is written as 



%-p(t) = i [H, p(t)\ - pT h p{t) - M r aP(*)> 
at in 



(3.6) 



where /1 = A 2 . In principle the form of the dissipation terms of Eq. ( 2.18| ) is derived in the condition where the system 
is coupled only to one reservoir . E ven when the two different reserv oirs a re contact with the system, the decoupled 
form of the dissipation term in (3.6) is valid in the order of A 2 . From Q2.1S| ) the damping terms and T R p(t) are 



T L p(t) = J {[ Xl ,R L p(t)] + [ Xl ,R h p(t)] j } and T R p{t) = J {[x N ,R R p(t)] + [x N . 
respectively. Here operators i?L and R R are defined through the matrix elements 



(l\R h \m) = 
(l\R R \m) = 



E,-E„ 



T (Ei—E m \ T I Ei-E„ 



^(Et-E,, 



1 



(Z|xi|to), 
(l\x N \m) 



(3.7) 

(3.8a) 
(3.8b) 



and Ei denotes the energy eigenvalue for state 7l and Ir are the spectral density of the left and the right reservoir, 
respectively. 

To solve this equation, we need to express the operators i?L and R R in terms of a,k and a\. The operators x\ and 
xjv are written as 



xi 



E 



sin k 



2{N + l)mwo ^ ^/ s in{k/2) 



(a k + 4) 



(3.9) 
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and 



x N 



' h ^ sm (Nk) r f 
2(N + l)mu \ v /sin(fc/2) V°* + ° fe 



Let denote the eigenstate for the number operator a\.a k with the eigenvalue n^, namely a\,a k \n k ) = n k 
eigenstates for the system Hamiltonian ( |3.5[ ) are given by the direct product of number states |n/-) as = 
whose energy eigenvalue is E({n k }) = Y^,k \ nk + i) 

The matrix elements of i?L are given in terms of the eigenstates as, 



(3.10) 
n k ). The 



({ny}\R h \{my}) 



2(N + l)mcj Y v/sinWS) 



E 



sin fc 



J L / g({n fc /})-g({m fc ,» 



-/l " 



Mi 



exp{/3 L [£({n fc ,}) - £({m fc >})]} - 1 



Now we note that 



[({nfc'IMWfc'}) + {{ n k'}W k \{my}) 



({ny}\a k \{my}) ^ only if ny = my - 8y ,k for Vfc', 
{{ n k'}\a\\{my}) ^ only if ny = my + S k ^ k for Vfc', 



and /(w) = for u < 0. Then Eq. (3.11) is transformed into 

h(u)k) sinfc 



({ny}\R h \{my}) = 



2(N+ l)mw Y VshiO/ 2 ) 



X 



1 



Thus the operator i?L can be represented as 

/ ?i ^ sm k 

8(N + l)muj \ 1 /sin(fc/2) sinh (fchuk/2) 



{{ny}\a\\{my}) 



Rl — 



In the same manner, the operator i?R is represented as 

sin(iVfc) 



Rr = 



E 



8(N + l)mu ^ v / sin(fc/2) sinh (j3 R hw k /2) 



(3.11) 

(3.12) 
(3.13) 



(3.14) 



(3.15) 



(3.16) 



C. Moments in the Stationary State 



As will be shown later, to evaluate mean kinetic energy of a particle and energy flux in the stationary state, we 
have only to calculate the second moments 



and 



{a k ay) = Tr(a k ayp st ) 



{a\ay) = Tr(a\.ayp s t), 



(3.17) 



(3.18) 



where p st denotes the stationary solution of Eq. fl3.6p . First, we consider Eq. (3.17). Because the left hand side of Eq. 
( |3.6|) vanishes in the stationary state, we obtain 



1 



7T/1 



Tr (a k a k > [H,p u 
in a 

7T/i 

~ IT 



Tr (a k ay [x u Rhpst}) + Tr (a k ay [x u i? L /O st ] t 
Tr (a k ay [x N ,R R p]) + Tr [a k ay [x N , i? R p] t 



= 0. 



(3.19) 
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This equation is rewritten as follows after tedious but straightforward calculations, 



i(uk + uj k >)(a k a k ,) + 
sin k' 



v /sin(fc'/2) 
sin k 



TT^ \ - sin fci J L (u kl ) 

4(7V + l)mm ^ VsR^72! sinh(/3 L ^ fel /2) 

=/L ^ l/2((afcafei) _ (4 i a fe )) + e-^^/ 2 ({a fc 4 i ) 



+ 



v /sin(A;/2) 

7T/U 



[e^/ 2 (K,a fcl > - (4 i a^)) + e-^/ 2 «a^4 i ) 
y-, sin(Affci) /R^feJ 



(a fcl a fc )) 
{a kl a k >}) 



4:(N+l)muj y/smih/2) sinh(/3 R ^w fcl /2) 



+ 



J sin(TVfc') 
\ Vsin(fc'/2) 
sin(iVfc) 
A/sin(A;/2) 



feR ^ /2 ((a fe a fel )-(4 i a fe ))+e-^ 1 / 2 ((a fe 4 i ) 
e fe««*x/2(( afc<0fci ) _ (4 i0fc ,)) +e-^^ 1 / 2 «a fc ,4 i ) ■ 



(a kl a k )) 
{a kl a k >)) 



In the same way, Eq. ( |3.18| ) is also transformed into 



i (w fe / - uj k ) {a\a k ,) + 



nil 



E 



sin fci 



4 (wfci ) 



A{N + l)mwo ^ v/shi(fcI72) sinh(/3 L 7iw fel /2) 



sin k 



+ 



{ v / sin(fc/2) 

sin fc' 
^sin(/c'/2) 



e' 3L ^ 1 / 2 «4 i a fe ,) - (a k ,a kl )) + e-^ h ^/ 2 ((a kl a k ,) - (a^)) 
[e^ R ^/ 2 ((4a fcl ) - (4 x 4)) + e-* R ^ 1 / 2 ((44 i ) - (a fcl 4» 



+ 



7T/i 



E 



sin(iVfci 



^R.( w fci) 



4(iV + l)mw ^ >/sin(A;i/2) sinh(^ R ftw fcl /2) 



J sin(TVfc) 
\ Vsin(fc/2) 

sin(iVfc') 
Vsin(A:'/2) 



e^ n ^ /2 ((ala k ,} - (a k ,a kl )) + e-^ h ^/ 2 ((a kl a k > 

^ 1 / 2 ((4a fcl ) - (4 1 4))+e- feR ^ /2 ((44 1 > - 



- (afe' a L>) 
Ki4>) 



0. 



(3.20) 



(3.21) 



D. Total Energy Est 

We will solve these equations by perturbation. Expanding (a k a k >) and (a k a k >) with respect to /i, 

(a*afc') = (afcafc')o + At(afcafe')i + M 2 ( a fc a fc'>2 + ■ 

(44-) = (44')o + m<44')i + /* 2 <44')a + • 

(4 a fc') = (4 a fe')o + ^( a t a k')i + H 2 ( a k a k>)2 + ■ 
{a k a\,) = (o fc 4'>o + ^(ofc4')i + M 2 (a/ta[./)2 + ■ 
we consider the relation of each order of /x. Using the commutation relations 

(a k a k ') n — (a k >a k )n, 

(44' )« = (4'4)«> 



(3.22) 
(3.23) 
(3.24) 
(3.25) 



(3.26) 
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we obtain the following relations at the zeroth order, 

(uj k + u] k ,)(a k a k ')o = and (u k > - w fe )(a[a fe /)o = 0. 
Accordingly we have for all k and k! 

(a k a k >) = 0, 

and for k =/= k' 

(%ah')o = 0. 

Putting k = k' in the first order equation of /z, we have 

(e^^ /2 (4« fc )o - e-^(a k a{) ) 
smh(/3 L ?iWfe/2) V K / 

smh(/3 R ^w fc /2) V 7 

Since sin 2 k = sin 2 (Nk) ^ 0, Eq. ( fT3C| ) leads to 

/ t \ 1 h(uJk) In(uk) 

I L (u k ) + I R {uj k ) 

Here the energy of the system is, 

flLO k 



E st =Tr{Hp st )=J2 



J L (w fe ) Ir(u))c) 



(3.27) 
(3.28) 
(3.29) 



(3.30) 



(3.31) 



(3.32) 



In particular, when Ii,(ui k ) = lR,(u)k), we find that E st is the arithmetic mean between equilibrium energy at inverse 
temperature /3l and at /3r regardless of the types of the spectral density, i.e. 



E R t — — 



Tr (He-^ H ) Tr (H e -^ H ) 



Tr e~fc H 



Tr e~^ H 



(3.33) 



E. Kinetic Energy of a Particle and Energy Flux 



Here we compute mean kinetic energy of each particle and energy flux up to the first order with respect to fi. The 
mean kinetic energy of the nth particle, e n is defined by 



Pn 

2m 



2 m 



which is expressed in terms of the creation and annihilation operators as 



N - 



- \J sin(fc/2) sin(fc'/2) sin(fcn) sin(fc'r 



k,k' 



x Ua k a k >) - {a k a\,) - (a\a k >) + {a\a\,) 
Substituting the results obtained in the last subsection into the above equation, we have 



N + 

hu> 



- ^2 sin(fc/2) sin 2 (kn) ^2(a^a fc ) + 



- sin(/c/2) sin 2 (kn) 



h(uk) 



-TR(wfc) 



h(u k ) + lR.(uk) \e^ hu " - 1 e ^ hu "> - 1 



(3.34) 



(3.35) 



(3.36) 







In the classical limit, Eq. ( 3.36| ) becomes 



2T L r . 2 h(uk) 2T R C „ . 2n \ lR{u k ) 
2e n = / dk sm (kn)-— — r — —. — - H / dk sin (kn)-—. — r — -. — -, (3.37) 

7T Jo ' h(w k ) + I K {w k ) 7T J I L {uj k ) + I R (u>k) 

where Tl = /3r 1 and Tr = (3^ 1 are the temperature values of the left reservoir and the right reservoir, respectively. 
In the classical limit, 2e„ can be interpreted as the temperature at site n. Especially when both the reservoirs are 
of the same type, namely, Ih(uJk) = iR^k), we have 2e n = Tl 1, Tr regardless of the types of the spectral density. 
This means completely flat temperature profile, which was originally found by RLL when the both reservoirs are of 



the Ohmic type ||. On the other hand, Eq. (3.37) shows that the temperature profile in integrable systems can 
be easily changed by controlling the combination of the types of spectral density of reservoirs, so that in the case of 
Ih(^k) 7^ I~R,(uk), the internal temperature deviates from Tl + Tr in the classical limit. 



We numerically estimate (3.36) to investigate the general feature of temperature profile for various combinations 
of spectral densities. We present typical temperature profiles in Fig.l. In Fig. 1(a) we take the sub-Ohmic type 
reservoir Tl = w °' 5 for left side, and the super-Ohmic one Ir = w 15 for the right side. In Fig. 1(b), the converse case, 
namely, Tl = w 1 ' 5 and Jr = oj 5 is considered. In both the cases, parameters are set to m = h = ujq = 1.0, and 
Tl = 200.0, Tr = 50.0. As is known from these figures, temperature deviates from the internal temperature value in 
the same direction in the vicinity of both the ends. As the result the deviated temperature values become close to 
the temperature of the reservoir whose spectral density has larger power. We numerically confirmed this dependence 
for many sets of spectral densities (-/l,Zr) and temperatures (Tl,Tr) including low temperatures. 

Energy flux is defined via the equation of continuity. From the master equation ( |3.6| ) , time derivative for the energy 
of the system satisfies 

J^Tr {Hp(t)) = -Tr (pHT^t)) - Tr {fiHT K p(t)) . (3.38) 

The first term in the right-hand side is regarded as incoming energy flux from the left reservoir and the second term 
incoming energy flux from the right reservoir. We call the former Jl and the latter Jr. In the stationary state, of 
course Jf* + JjJ = must hold. We can calculate J£* as follows, 



(N + l)m Y Z L (w fc ) + IrK) 
If N 3> 1, we can replace the summation by an integral and obtain 

™t hf* r Z L (w fe )/R(w fe ) 2,7 1 ! 



(3.39) 



JL =— / r / \ V / \ sin J k t= - -s-s 7 dk. (3.40) 

In the classical limit (% — > 0), J£* goes to 

J£ 4 = M C(T L - Tr), (3.41) 

where 

c _ 1 f sin 2 fc I L (2^ sin(fc/2))J R (2u; sin(fc/2)) ^ 

mujo Jq sin(fc/2) / L (2w sin(fc/2)) + Ir(2w sin(fc/2)) ' [ ' ' 

Thus in the classical limit, energy flux is proportional to the temperature difference and independent of the system 
size regardless of the types and the combinations of the spectral densities of reservoirs. 



IV. FINITE COUPLING 



The master equation (2.18) is justified only in O(fi). However when we study the model with a finite coupling 
constant, the quantitative effect of a finite coupling inevitably deviates from those of the original model. Nevertheless, 
time evolution of the reduced density matrix has been regarded as describing a variety of relaxation processes, and it 
successfully explained a variety of interesting phenomena in real systems ]l7|,|l9| • This shows that the master equation 
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can well approximate the dynamics in real dissipative system at least qualitatively. In some case the time evolution 
of reduced density matrix with a finite coupling reproduces quantitatively correct results even for long time [fl0| . 

Thus we investigate here the effect of finite coupling which is a small but finite value, and discuss the behavior of 
temperature profile qualitatively. In this section, we confine ourselves to the case of the same spectral density at both 
ends, namely 

I L =I R = I. (4.1) 



A. Temperature Profile 



We evaluate contributions from higher-order terms and find deviations from the flat temperatur e pro file ne ar th e 
ends of the chain. We first calculate the first-order coefficients. From the first-order equations in (|3.20| ) and ( 3.21 ), 
we have for all k and k' 



(akO-k')i 



iir [sin k sin k' — s'm(Nk) sm(Nk')] 



and for k =/= k' 



(4 a fc')l = 



4(N + l)muo(uJk + uJ k ')y/sm(k/2) sm(fc'/2) 

x [I(u) k ) (n L (w fe ) - n R (u> k )) + i(cJfe') (n L (wfc') - n R (u k '))] ; 

in [sin k sin k' — s'm(Nk) sin(iVfc')] 



(4.2) 



4(N + l)muJo{uJk - uJ k >)y / sm(k/2) sin(fc'/2) 
x [I(wfe) (n L (wfe) - n R (u> k )) + I(oJk>) (n h (uj k >) - n R (uj k >))} , 
where ni,(u>) is the Bose-Einstein distribution functions at an inverse temperature /3l 

1 



and n R (uj) is that at (3 R 



n L (w) 



n R (uj) 



1 



(4.3) 
(4.4) 

(4.5) 



Equations ( 3.20| ) and (3.21) imply that the first-order coefficients must be pure imaginary. On the other hand, (a k a k )i 
must be real at the same time. Thus, we have 



{a\a k )i = 0. 



(4.6) 

From Eqs. ( 3.20|) and ( 3.21 ), it turns out that if n > 1 the (n + l)st-order terms are computed via the following 
equations from the nth order terms; For all k and kl , 



(a k a k ') n 



+i 



2(JV + l)mwoK + ^) J2 T M 

k i 

sinfci sin kl + sin(TVfci) sin(iYfc') 
x /sin(fci/2)sin(fc'/2) 

sin ki sin k + sm(Nki) sin(iVfc) / t 

- ( \a k ia kl ) n - {a' ki a k <), 



((a k a kl ) n - (a^afe), 



+ 



If k ^ k', 



{a\a k ') n +i = 



A/sin(A;i/2)sin(fc/2) 



(4.7) 



2(N + l)mLU (Lu k < - uj k ) 



sin ki sin kl + sin(iVfci ) sin(iVfc') 
v / sin(fci/2)sin(fc'/2) 

sinfci sin A; + sin(A r fc 1 ) sin(iVfc) 
v /sin(/ci/2)sin(fc/2) 



a k a kl ) n - {aka kl ) n 
( a kl a k')n - (a k >a kl ) 



(4.8) 
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and (alcik)n+l is computed through the other coefficients of the same order as 

(<4fflk)n+l = 2 ((°fcafc)n+i + (a k a k )* l+1 ) 

Jsm(k/2) ^ sin k' sin fc + sinfiVA:') sin(TVfc) , 

/ ■ frfm T M 



4/(cj fe )sin fc^ #fc 



v /sin(fc'/2) 

f(<4>afc)n+i + (a{a k ') n +i ~ (afeafe')n+i - (afc«fc 



)n+l) • 



(4.9) 



Becau se the above equations contain the spectral density, we have to specify its functional form. As has given in 
(2.13), we employ the following form for the spectral density 



J(W) = I LO a . 



(4.10) 



For a — 1.0 an d a — 1.5, we have computed mean kinetic energy of the nth particle e n up to the 20th order, where 
quantities ( |3.22 )-( 3.25D seem to converge. For each a, the following sets of temperatures of for the reservoirs are 
chosen: (a) T L = 200.0 and T R = 50.0, (b) T L = 10.0 and T R = 0.1 and (c) T L = 0.1 and T R = 0.02. These choices 
of parameters are the same as used in Q by ZT. Other parameters are commonly set as m = loq = Ti = 1.0, /i = 0.1, 
and Iq — 1/V. In the equilibrium state at inverse temperature j3, s n is given by 



MP) 



hujo 
ivT 



N 

1=1 



nl 



ulr 



sin 



sm 



2(JV+1) N+l 



■ coth 



(3fiujQ sin 



2(JV+1) 



(4.11) 



and thus the local temperatures T n is defined by the above function, i.e. T n = l/(f>~ 1 (s n ). 

In fig. 2 and 3, {T n } are plotted for a = 1.0 and 1.5, respectively. All the figures show that higher-order contributions 
are small except near the ends of the chain. In other words, the bulk behavior is unchanged, where the temperature 
profile near the ends of the chain exhibits various dependencies on the details of the parameters (Xl, Tr, a). When 
the reservoirs are Ohmic, the temperature profile near the ends are similar to those obtained by ZT with the quantum 
Langevin approach. 

When the reservoirs are Ohmic and temperature is high (Fig. 2(a)), temperature drops near the left end which 
contacts with the hotter reservoir and rises near the other end contacting with the colder reservoir. This is the same 
paradoxical behavior as found by RLL and also observed by ZT. Such behavior disappear when the reservoirs are 
super-Ohmic (Fig. 3(a)). The second particles from the ends shows monotonic temperature variation. 

Figure 2 (c) and Fig. 3 (c) exhibit temperature profile when the temperatures are low where quantum effects are 
important. These two figures almost coincide. In both cases, temperatures near the ends are high which should be due 
to the quantum fluctuations. We may say that differences in the spectral densities does not affect the temperature 
profile at low temperatures. In the medium temperature cases, Fig. 2 (b) and Fig. 3 (b), mixed behavior of the 
classical and quantum features are observed. (See also the figures in reference Q). 

For XL = 200.0 and Tr = 50.0, temperature deviations of the particle 2 and the particle (N — 1) from the mean 
internal temperature are plotted in Fig. 4 for various a. There we find that the peculiarity, i.e. inversion of temperature 
near the ends, is observed in the sub-Ohmic and Ohmic cases, while it disappears when a > a c {— 1-04). 



B. Fokker-Planck Equation in the Classical Limit 



In the previous subsection, the temperature profiles is found to depend on values of a. In particular, the peculiarity 
found by RLL || disappears in the super-Ohmic regime. Since the differences are seen at high temperatures, some 
characteristics depending to the values a must appear in the Fokker-Planck equation obtained from the master 
equation in the classical limit. Actually, we will find that the diffusion term in the Fokker-Planck equation takes a 
different form from that derived from the Langevin equation except the Ohmic case. In order to study the difference 
in relaxation at the contacting point, we investigate the Fokker-Planck equation for a system with a single reservoir. 

When a heat reservoir is attached to the left end of the chain, the classical Langevin equations for canonical variables 
x n (t) and p n (t), n — 1, 2, N are 



= {p n ,H} ~ 6 n}1 V— + 6 n ,l£(t) 

at m 



(4.12) 
(4.13) 
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where {•, •} means the Poisson bracket. The correlation function of the Gaussian white random force is connected 
with the damping constant v and the temperature at the first particle (3 via the fluctuation-dissipation theorem as 



mm) = jKt-t'). 

As is well known, the Langevin equations are equivalent to the Fokker-Planck equation 



dP(t) 
dt 



= {H,P(t)}+u 



d (pi 



dpi \ m 



dpi 



(4.14) 



(4.15) 



where P(t) is the distribution function on the phase space. 

W e now turn to our master equation. Inserting the representation for operator i?L ( 3.15 ) into the master equation 
([2.6D, we have 



7T/Z 



sin k 



dp(t) = 1 m (t)] __ 

dt ih [ ,P{>1 ^8(7V + rjmtkJo ^sin(fc/2) sinh(/3?iw fc /2) 



(4.16) 



{ [ Xl , (e^»' 2 a k + e-^/ 2 al) p(t)] - \x u p{t) (e^ 2 a\ + e~^ 2 a k )\ } , 



where we omitted suffix L. Expressing the creation and annihilation operators by the position and momentum oper- 
ators, 



dp® 
dt 



\[H, P (t)} 



1 (uj k ) sin ksin(kn) jcoth (^fl^j [ Xji) 



Pit)]] 



2mwo sin(fc/2) 



(Pn[xi,p(t)) + [X 1 ,p{t)]p n + 2[xi, p n ]p(t)) 



(4.17) 



In the classical limit, the density matrix p(t) is replaced by the distribution function P(t) Therefore, Eq. ( 4.17 ) is 
transformed into 



dP(t) 
~dT 



N 



{#,P(t)} + £c„ 



d ( p 



dpi \ m 



dpr, 



P(t), 



(4.18) 



where 



2p 
^0 Jo 



I(2uj sin(fe/2)) cos(fc/2) sm(kn)dk, 



(4.19) 



The time-evolution equation for the covariance matrix derived from Eq. (|4.18j ) is also confirmed to agree with the 
classical limit of the corresponding quantum equation. 

When the reservoir is Ohmic, namely I(uj) = Iqlu, the coefficients C n are evaluated as 



C n = npIoSn,!, 



(4.20) 



and Eq. ( 4.18 ) agre es wi th the Fokker-Planck equation derived from the Langevin equation ( 4.15 ). In this case, the 
two-point function ( 2.1Cf ) tends to the delta function in the classical limit 



Bm*(*) = ^(t). 



(4.21) 



Therefore, we find t hat t he correlation function of the noise is white in the Ohmic case, which is consistent with the 
Langevin equation ( |L13|) . 

If the reservoir is sub-Ohmic or super-Ohmic, however, C n does not vanish for n > 2. Figure 5 shows C n as a 
function of n for various values of a. The sign of C n (n > 2) is positive in the sub-Ohmic regime and negative in the 
super-Ohmic regime. The difference in temperature profiles discussed in the previous subsection should be explained 
by this a-dependence of the coefficients C n . 
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V. SUMMARY AND DISCUSSION 



We investigated the effect of the types of reservoirs on the thermal conduction in the harmonic chain. We derived 
the master equation for general many body system in contact with phonon reservoirs. In many body system, the 
dissipation term is different from one of one-particle system due to the noncommutability of many body interaction, 
s o tha t the dissip ation term has rather complicated form. However we have the explicit form for the dissipation term 
( 2.16 ) and ( 2.18 ). We used it as the basic equation to study behavior of the system. The equation generally satisfies 
the necessary condition for the master equation that the canonical distribution must be a stationary solution when 
the reservoirs are at the same temperature. 

In Sec. 3, we have applied the master equation to energy transport in the quantum harmonic chain. We attached a 
phonon reservoir at one end and another at the other end. At weak coupling limit (A — > 0), we obtained explicit form 
of internal energy and energy flux. We rigorously proved that when the spectral densities of the reservoirs are of same 
type, the total energy of the system takes the arithmetic mean of the equilibrium energies at Tl and Tr regardless 
of the type of the spectral density. This result leads to the classical temperature Tl "^ Tr which is originally found by 
RLL using the Ohmic type of reservoir. On the other hand, when the types of spectral densities are different, the 
internal temperature is a function of the both densities, so that the temperature does not converge to Tl "^ Tr in the 
classical limit. The difference of spectral densities induces the deviation of temperatures around the both edges from 
the internal value. The temperature in the vicinity of both ends become close to the temperature of the reservoir 
whose spectral density has larger power. We numerically confirmed that this feature is general when the reservoirs 
are of different types. 

We numerically investigated the effect of finite coupling. We considered only the case of the same spectral densities 
of the reservoirs. Finite coupling contributes to the temperature profile only near the ends of the chain and bulk 
behavior is the same as that of weak coupling limit. We found that the profile near the ends depends on the spectral 
density for the reservoirs. When the reservoir temperature is sufficiently low where quantum fluctuations are dominant, 
temperature growth near both the ends was observed in every case. When the reservoir temperatures are high enough 
and the reservoir is sub-Ohmic or Ohmic, the same peculiar behavior, i.e., nonmonotonic change of the temperature, 
is observed as found in p[. However, in the case of super-Ohmic reservoir, the peculiarity disappears. 

In order to understand the dependence on the spectral density, we derived Fokker-Planck equations from the 
master equation in the classical limit. If the reservoir is Ohmic the Fokker-Planck equation agrees with the standard 
one derived from the Langevin equation. When the reservoir is non-Ohmic, however, there appears difference in the 
diffusion term , i.e., form of second derivative. The coefficients of the diffusion terms were calculated from the spectral 
density.This difference causes different temperature profiles near the ends of the chain. 

We expect that the master equation derived here can be used for other systems such as spin systems for which the 
Langevin approach is practically difficult. In the case of the harmonic chain, operator R was written in a simple form 
by using some system operators. Thus we were able to analyze the master equation systematically. This can be done 
because the harmonic chain is integrable. Thus, similar procedure can be developed for other integrable systems, e.g. 
the XY model §§. 

In this paper we have confined ourselves to the integrable system. However the master equation derived here is 
generally applicable to any system because the matrix element of R is explicitly given. Thus it would be an interesting 
future problem to study the thermal conductivity in nonintegrable system where the Fourier heat law is realized . 
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FIG. 1. Temperature profile along the chain for T L = 200.0, T R = 50.0: (a) J L (w) = c<J ' 5 and = a; 1 ' 5 .; (b) J L (w) = w 1 " 8 

and 7r(w) = u> ' 8 . The system size is N = 150. 



FIG. 2. Temperature profile along the chain for a = 1.0: (a) T L = 200.0, Tr = 50.0; (b) T L = 10.0, Tr = 0.1; (c) 
T L = 0.1, T R = 0.02. The system size is N = 150. 



FIG. 3. Temperature profile along the chain for a = 1.5: (a) T L = 200.0, Tr = 50.0; (b) T L = 10.0, Tr = 0.1; (c) 
T L = 0.1, Tr = 0.02. The system size is N = 150. 



FIG. 4. Deviations of the temperature at particle 2 and particle (N — 1) from the mean internal temperature T av . The 
temperatures of the reservoirs are Tl = 200.0 and Tr = 50.0. Thus, T av = 125.0. 



FIG. 5. Coefficients C n as a function of n for various values of a. 
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